#ifndef VECTOR4_HPP_
#define VECTOR4_HPP_

#include <cmath>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <string>

#include "Vector3.hpp"

class LorentzRotation;

class Vector4 {

private:

   Vector3 fP;  // 3 vector component
   double fE;  // time or energy of (x,y,z,t) or (px,py,pz,e)

public:

   enum { kX=0, kY=1, kZ=2, kT=3, kNUM_COORDINATES=4, kSIZE=kNUM_COORDINATES };
   // Safe indexing of the coordinates when using with matrices, arrays, etc.

   Vector4(double x = 0.0, double y = 0.0,
                 double z = 0.0, double t = 0.0);
   // Constructor giving the components x, y, z, t.

   Vector4(const double * carray);
   Vector4(const float * carray);
   // Constructor from an array, not checked!

   Vector4(const Vector3 & vector3, double t);
   // Constructor giving a 3-Vector and a time component.

   Vector4(const Vector4 & lorentzvector);
   // Copy constructor.

   virtual ~Vector4();
   // Destructor

   // inline operator Vector3 () const;
   // inline operator Vector3 & ();
   // Conversion (cast) to Vector3.

   inline double X() const;
   inline double Y() const;
   inline double Z() const;
   inline double T() const;
   // Get position and time.

   inline void SetX(double a);
   inline void SetY(double a);
   inline void SetZ(double a);
   inline void SetT(double a);
   // Set position and time.

   inline double Px() const;
   inline double Py() const;
   inline double Pz() const;
   inline double P()  const;
   inline double E()  const;
   inline double Energy() const;
   // Get momentum and energy.

   inline void SetPx(double a);
   inline void SetPy(double a);
   inline void SetPz(double a);
   inline void SetE(double a);
   // Set momentum and energy.

   inline Vector3 Vect() const ;
   // Get spatial component.

   inline void SetVect(const Vector3 & vect3);
   // Set spatial component.

   inline double Theta() const;
   inline double CosTheta() const;
   inline double Phi() const; //returns phi from -pi to pi
   inline double Rho() const;
   // Get spatial vector components in spherical coordinate system.

   inline void SetTheta(double theta);
   inline void SetPhi(double phi);
   inline void SetRho(double rho);
   // Set spatial vector components in spherical coordinate system.

   inline void SetPxPyPzE(double px, double py, double pz, double e);
   inline void SetXYZT(double  x, double  y, double  z, double t);
   inline void SetXYZM(double  x, double  y, double  z, double m);
   inline void SetPtEtaPhiM(double pt, double eta, double phi, double m);
   inline void SetPtEtaPhiE(double pt, double eta, double phi, double e);
   // Setters to provide the functionality (but a more meanigful name) of
   // the previous version eg SetV4... PsetV4...

   inline void GetXYZT(double *carray) const;
   inline void GetXYZT(float *carray) const;
   // Getters into an arry
   // no checking!

   double operator () (int i) const;
   inline double operator [] (int i) const;
   // Get components by index.

   double & operator () (int i);
   inline double & operator [] (int i);
   // Set components by index.

   inline Vector4 & operator = (const Vector4 &);
   // Assignment.

   inline Vector4   operator +  (const Vector4 &) const;
   inline Vector4 & operator += (const Vector4 &);
   // Additions.

   inline Vector4   operator -  (const Vector4 &) const;
   inline Vector4 & operator -= (const Vector4 &);
   // Subtractions.

   inline Vector4 operator - () const;
   // Unary minus.

   inline Vector4 operator * (double a) const;
   inline Vector4 & operator *= (double a);
   // Scaling with real numbers.

   inline bool operator == (const Vector4 &) const;
   inline bool operator != (const Vector4 &) const;
   // Comparisons.

   inline double Perp2() const;
   // Transverse component of the spatial vector squared.

   inline double Pt() const;
   inline double Perp() const;
   // Transverse component of the spatial vector (R in cylindrical system).

   inline void SetPerp(double);
   // Set the transverse component of the spatial vector.

   inline double Perp2(const Vector3 & v) const;
   // Transverse component of the spatial vector w.r.t. given axis squared.

   inline double Pt(const Vector3 & v) const;
   inline double Perp(const Vector3 & v) const;
   // Transverse component of the spatial vector w.r.t. given axis.

   inline double Et2() const;
   // Transverse energy squared.

   inline double Et() const;
   // Transverse energy.

   inline double Et2(const Vector3 &) const;
   // Transverse energy w.r.t. given axis squared.

   inline double Et(const Vector3 &) const;
   // Transverse energy w.r.t. given axis.

   inline double DeltaPhi(const Vector4 &) const;
   inline double DeltaR(const Vector4 &) const;
   inline double DrEtaPhi(const Vector4 &) const;
   inline Vector2 EtaPhiVector();

   inline double Angle(const Vector3 & v) const;
   // Angle wrt. another vector.

   inline double Mag2() const;
   inline double M2() const;
   // Invariant mass squared.

   inline double Mag() const;
   inline double M() const;
   // Invariant mass. If mag2() is negative then -sqrt(-mag2()) is returned.

   inline double Mt2() const;
   // Transverse mass squared.

   inline double Mt() const;
   // Transverse mass.

   inline double Beta() const;
   inline double Gamma() const;

   inline double Dot(const Vector4 &) const;
   inline double operator * (const Vector4 &) const;
   // Scalar product.

   inline void SetVectMag(const Vector3 & spatial, double magnitude);
   inline void SetVectM(const Vector3 & spatial, double mass);
   // Copy spatial coordinates, and set energy = sqrt(mass^2 + spatial^2)

   inline double Plus() const;
   inline double Minus() const;
   // Returns t +/- z.
   // Related to the positive/negative light-cone component,
   // which some define this way and others define as (t +/- z)/sqrt(2)

   inline Vector3 BoostVector() const ;
   // Returns the spatial components divided by the time component.

   void Boost(double, double, double);
   inline void Boost(const Vector3 &);
   // Lorentz boost.

   double Rapidity() const;
   // Returns the rapidity, i.e. 0.5*ln((E+pz)/(E-pz))

   inline double Eta() const;
   inline double PseudoRapidity() const;
   // Returns the pseudo-rapidity, i.e. -ln(tan(theta/2))

   inline void RotateX(double angle);
   // Rotate the spatial component around the x-axis.

   inline void RotateY(double angle);
   // Rotate the spatial component around the y-axis.

   inline void RotateZ(double angle);
   // Rotate the spatial component around the z-axis.

   inline void RotateUz(Vector3 & newUzVector);
   // Rotates the reference frame from Uz to newUz (unit vector).

   inline void Rotate(double, const Vector3 &);
   // Rotate the spatial component around specified axis.

   inline Vector4 & operator *= (const Rotation &);
   inline Vector4 & Transform(const Rotation &);
   // Transformation with HepRotation.

   Vector4 & operator *= (const LorentzRotation &);
   Vector4 & Transform(const LorentzRotation &);
   // Transformation with HepLorenzRotation.

  std::string ToString();

};


class LorentzRotation {

public:

#ifndef SWIG
class LorentzRotationRow {
public:
   inline LorentzRotationRow(const LorentzRotation &, int);
   inline double operator [] (int) const;
private:
   const LorentzRotation * fRR;
   int fII;
};
#endif
   // Helper class for implemention of C-style subscripting r[i][j]

   LorentzRotation();
   // Default constructor. Gives a unit matrix.

   LorentzRotation(const Rotation &);
   // Constructor for 3d rotations.

   LorentzRotation(const LorentzRotation &);
   // Copy constructor.

   LorentzRotation(double, double, double);
   LorentzRotation(const Vector3 &);
   // Constructors giving a Lorenz-boost.

   inline double XX() const;
   inline double XY() const;
   inline double XZ() const;
   inline double XT() const;
   inline double YX() const;
   inline double YY() const;
   inline double YZ() const;
   inline double YT() const;
   inline double ZX() const;
   inline double ZY() const;
   inline double ZZ() const;
   inline double ZT() const;
   inline double TX() const;
   inline double TY() const;
   inline double TZ() const;
   inline double TT() const;
   // Elements of the matrix.

   inline LorentzRotationRow operator [] (int) const;
   // Returns object of the helper class for C-style subscripting r[i][j]


   double operator () (int, int) const;
   // Fortran-style subscriptimg: returns (i,j) element of the matrix.


   inline LorentzRotation & operator = (const LorentzRotation &);
   inline LorentzRotation & operator = (const Rotation &);
   // Assignment.

   inline bool operator == (const LorentzRotation &) const;
   inline bool operator != (const LorentzRotation &) const;
   // Comparisons.

   inline bool IsIdentity() const;
   // Returns true if the Identity matrix.

   inline Vector4 VectorMultiplication(const Vector4&) const;
   inline Vector4 operator * (const Vector4 &) const;
   // Multiplication with a Lorentz vector.

   LorentzRotation MatrixMultiplication(const LorentzRotation &) const;
   inline LorentzRotation operator * (const LorentzRotation &) const;
   inline LorentzRotation & operator *= (const LorentzRotation &);
   inline LorentzRotation & Transform(const LorentzRotation &);
   inline LorentzRotation & Transform(const Rotation &);
   // Matrix multiplication.
   // Note: a *= b; <=> a = a * b; while a.Transform(b); <=> a = b * a;

   inline LorentzRotation Inverse() const;
   // Return the inverse.

   inline LorentzRotation & Invert();
   // Inverts the LorentzRotation matrix.

   inline LorentzRotation & Boost(double, double, double);
   inline LorentzRotation & Boost(const Vector3 &);
   // Lorenz boost.

   inline LorentzRotation & RotateX(double);
   // Rotation around x-axis.

   inline LorentzRotation & RotateY(double);
   // Rotation around y-axis.

   inline LorentzRotation & RotateZ(double);
   // Rotation around z-axis.

   inline LorentzRotation & Rotate(double, const Vector3 &);
   inline LorentzRotation & Rotate(double, const Vector3 *);
   // Rotation around specified vector.

  std::string ToString();

protected:

   double fxx, fxy, fxz, fxt,
            fyx, fyy, fyz, fyt,
            fzx, fzy, fzz, fzt,
            ftx, fty, ftz, ftt;
   // The matrix elements.

   void SetBoost(double, double, double);
   // Set elements according to a boost vector.

   LorentzRotation(double, double, double, double,
                    double, double, double, double,
                    double, double, double, double,
                    double, double, double, double);
   // Protected constructor.

};



//inline Vector4 operator * (const Vector4 &, double a);
// moved to Vector4::operator * (double a)
inline Vector4 operator * (double a, const Vector4 &);
// Scaling LorentzVector with a real number


inline double Vector4::X() const { return fP.X(); }
inline double Vector4::Y() const { return fP.Y(); }
inline double Vector4::Z() const { return fP.Z(); }
inline double Vector4::T() const { return fE; }

inline void Vector4::SetX(double a) { fP.SetX(a); }
inline void Vector4::SetY(double a) { fP.SetY(a); }
inline void Vector4::SetZ(double a) { fP.SetZ(a); }
inline void Vector4::SetT(double a) { fE = a; }

inline double Vector4::Px() const { return X(); }
inline double Vector4::Py() const { return Y(); }
inline double Vector4::Pz() const { return Z(); }
inline double Vector4::P()  const { return fP.Mag(); }
inline double Vector4::E()  const { return T(); }
inline double Vector4::Energy()  const { return T(); }

inline void Vector4::SetPx(double a) { SetX(a); }
inline void Vector4::SetPy(double a) { SetY(a); }
inline void Vector4::SetPz(double a) { SetZ(a); }
inline void Vector4::SetE(double a)  { SetT(a); }

inline Vector3 Vector4::Vect() const { return fP; }

inline void Vector4::SetVect(const Vector3 &p) { fP = p; }

inline double Vector4::Phi() const {
   return fP.Phi();
}

inline double Vector4::Theta() const {
   return fP.Theta();
}

inline double Vector4::CosTheta() const {
   return fP.CosTheta();
}


inline double Vector4::Rho() const {
   return fP.Mag();
}

inline void Vector4::SetTheta(double th) {
   fP.SetTheta(th);
}

inline void Vector4::SetPhi(double phi) {
   fP.SetPhi(phi);
}

inline void Vector4::SetRho(double rho) {
   fP.SetMag(rho);
}

inline void Vector4::SetXYZT(double  x, double  y, double  z, double t) {
   fP.SetXYZ(x, y, z);
   SetT(t);
}

inline void Vector4::SetPxPyPzE(double px, double py, double pz, double e) {
   SetXYZT(px, py, pz, e);
}

inline void Vector4::SetXYZM(double  x, double  y, double  z, double m) {
   if ( m  >= 0 ) 
      SetXYZT( x, y, z, sqrt(x*x+y*y+z*z+m*m) );
   else {
     double mass = x*x+y*y+z*z-m*m;
     SetXYZT( x, y, z, (mass > 0) ? mass : 0);
   }
}

inline void Vector4::SetPtEtaPhiM(double pt, double eta, double phi, double m) {
   pt = fabs(pt);
   SetXYZM(pt*cos(phi), pt*sin(phi), pt*sinh(eta) ,m);
}

inline void Vector4::SetPtEtaPhiE(double pt, double eta, double phi, double e) {
   pt = fabs(pt);
   SetXYZT(pt*cos(phi), pt*sin(phi), pt*sinh(eta) ,e);
}

inline void Vector4::GetXYZT(double *carray) const {
   fP.GetXYZ(carray);
   carray[3] = fE;
}

inline void Vector4::GetXYZT(float *carray) const{
   fP.GetXYZ(carray);
   carray[3] = fE;
}

double & Vector4::operator [] (int i)       { return (*this)(i); }
double   Vector4::operator [] (int i) const { return (*this)(i); }

inline Vector4 &Vector4::operator = (const Vector4 & q) {
   fP = q.Vect();
   fE = q.T();
   return *this;
}

inline Vector4 Vector4::operator + (const Vector4 & q) const {
   return Vector4(fP+q.Vect(), fE+q.T());
}

inline Vector4 &Vector4::operator += (const Vector4 & q) {
   fP += q.Vect();
   fE += q.T();
   return *this;
}

inline Vector4 Vector4::operator - (const Vector4 & q) const {
   return Vector4(fP-q.Vect(), fE-q.T());
}

inline Vector4 &Vector4::operator -= (const Vector4 & q) {
   fP -= q.Vect();
   fE -= q.T();
   return *this;
}

inline Vector4 Vector4::operator - () const {
   return Vector4(-X(), -Y(), -Z(), -T());
}

inline Vector4& Vector4::operator *= (double a) {
   fP *= a;
   fE *= a;
   return *this;
}

inline Vector4 Vector4::operator * (double a) const {
   return Vector4(a*X(), a*Y(), a*Z(), a*T());
}

inline bool Vector4::operator == (const Vector4 & q) const {
   return (Vect() == q.Vect() && T() == q.T());
}

inline bool Vector4::operator != (const Vector4 & q) const {
   return (Vect() != q.Vect() || T() != q.T());
}

inline double Vector4::Perp2() const  { return fP.Perp2(); }

inline double Vector4::Perp()  const  { return fP.Perp(); }

inline double Vector4::Pt() const { return Perp(); }

inline void Vector4::SetPerp(double r) {
   fP.SetPerp(r);
}

inline double Vector4::Perp2(const Vector3 &v) const {
   return fP.Perp2(v);
}

inline double Vector4::Perp(const Vector3 &v) const {
   return fP.Perp(v);
}

inline double Vector4::Pt(const Vector3 &v) const {
   return Perp(v);
}

inline double Vector4::Et2() const {
   double pt2 = fP.Perp2();
   return pt2 == 0 ? 0 : E()*E() * pt2/(pt2+Z()*Z());
}

inline double Vector4::Et() const {
   double etet = Et2();
   return E() < 0.0 ? -sqrt(etet) : sqrt(etet);
}

inline double Vector4::Et2(const Vector3 & v) const {
   double pt2 = fP.Perp2(v);
   double pv = fP.Dot(v.Unit());
   return pt2 == 0 ? 0 : E()*E() * pt2/(pt2+pv*pv);
}

inline double Vector4::Et(const Vector3 & v) const {
   double etet = Et2(v);
   return E() < 0.0 ? -sqrt(etet) : sqrt(etet);
}

inline double Vector4::DeltaPhi(const Vector4 & v) const {
   return Vector2::Phi_mpi_pi(Phi()-v.Phi());
}

inline double Vector4::Eta() const {
   return PseudoRapidity();
}
inline double Vector4::DeltaR(const Vector4 & v) const {
   double deta = Eta()-v.Eta();
   double dphi = Vector2::Phi_mpi_pi(Phi()-v.Phi());
   return sqrt( deta*deta+dphi*dphi );
}

inline double Vector4::DrEtaPhi(const Vector4 & v) const{
   return DeltaR(v);
}

inline Vector2 Vector4::EtaPhiVector() {
   return Vector2 (Eta(),Phi());
}


inline double Vector4::Angle(const Vector3 &v) const {
   return fP.Angle(v);
}

inline double Vector4::Mag2() const {
   return T()*T() - fP.Mag2();
}

inline double Vector4::Mag() const {
   double mm = Mag2();
   return mm < 0.0 ? -sqrt(-mm) : sqrt(mm);
}

inline double Vector4::M2() const { return Mag2(); }
inline double Vector4::M() const { return Mag(); }

inline double Vector4::Mt2() const {
   return E()*E() - Z()*Z();
}

inline double Vector4::Mt() const {
   double mm = Mt2();
   return mm < 0.0 ? -sqrt(-mm) : sqrt(mm);
}

inline double Vector4::Beta() const {
   return fP.Mag() / fE;
}

inline double Vector4::Gamma() const {
   double b = Beta();
   return 1.0/sqrt(1- b*b);
}

inline void Vector4::SetVectMag(const Vector3 & spatial, double magnitude) {
   SetXYZM(spatial.X(), spatial.Y(), spatial.Z(), magnitude);
}

inline void Vector4::SetVectM(const Vector3 & spatial, double mass) {
   SetVectMag(spatial, mass);
}

inline double Vector4::Dot(const Vector4 & q) const {
   return T()*q.T() - Z()*q.Z() - Y()*q.Y() - X()*q.X();
}

inline double Vector4::operator * (const Vector4 & q) const {
   return Dot(q);
}

//Member functions Plus() and Minus() return the positive and negative
//light-cone components:
//
//  double pcone = v.Plus();
//  double mcone = v.Minus();
//
//CAVEAT: The values returned are T{+,-}Z. It is known that some authors
//find it easier to define these components as (T{+,-}Z)/sqrt(2). Thus
//check what definition is used in the physics you're working in and adapt
//your code accordingly.

inline double Vector4::Plus() const {
   return T() + Z();
}

inline double Vector4::Minus() const {
   return T() - Z();
}

inline Vector3 Vector4::BoostVector() const {
   return Vector3(X()/T(), Y()/T(), Z()/T());
}

inline void Vector4::Boost(const Vector3 & b) {
   Boost(b.X(), b.Y(), b.Z());
}

inline double Vector4::PseudoRapidity() const {
   return fP.PseudoRapidity();
}

inline void Vector4::RotateX(double angle) {
   fP.RotateX(angle);
}

inline void Vector4::RotateY(double angle) {
   fP.RotateY(angle);
}

inline void Vector4::RotateZ(double angle) {
   fP.RotateZ(angle);
}

inline void Vector4::RotateUz(Vector3 &newUzVector) {
   fP.RotateUz(newUzVector);
}

inline void Vector4::Rotate(double a, const Vector3 &v) {
   fP.Rotate(a,v);
}

inline Vector4 &Vector4::operator *= (const Rotation & m) {
   fP *= m;
   return *this;
}

inline Vector4 &Vector4::Transform(const Rotation & m) {
   fP.Transform(m);
   return *this;
}

inline Vector4 operator * (double a, const Vector4 & p) {
   return Vector4(a*p.X(), a*p.Y(), a*p.Z(), a*p.T());
}

Vector4::Vector4(double x, double y, double z, double t)
               : fP(x,y,z), fE(t) {}

Vector4::Vector4(const double * x0)
               : fP(x0), fE(x0[3]) {}

Vector4::Vector4(const float * x0)
               : fP(x0), fE(x0[3]) {}

Vector4::Vector4(const Vector3 & p, double e)
               : fP(p), fE(e) {}

Vector4::Vector4(const Vector4 & p) : fP(p.Vect()), fE(p.T()) {}

Vector4::~Vector4()  {}

double Vector4::operator () (int i) const
{
   //dereferencing operator const
   switch(i) {
      case kX:
      case kY:
      case kZ:
         return fP(i);
      case kT:
         return fE;
      default:
	std::cerr << "Vector4::operator()() bad index (" << i << "), returning 0" << std::endl;
   }
   return 0.;
}

double & Vector4::operator () (int i)
{
   //dereferencing operator
   switch(i) {
      case kX:
      case kY:
      case kZ:
         return fP(i);
      case kT:
         return fE;
      default:
	std::cerr << "Vector4::operator()() bad index (" << i << "), returning &fE" << std::endl;
   }
   return fE;
}

void Vector4::Boost(double bx, double by, double bz)
{
   //Boost this Lorentz vector
   double b2 = bx*bx + by*by + bz*bz;
   register double gamma = 1.0 / sqrt(1.0 - b2);
   register double bp = bx*X() + by*Y() + bz*Z();
   register double gamma2 = b2 > 0 ? (gamma - 1.0)/b2 : 0.0;

   SetX(X() + gamma2*bp*bx + gamma*bx*T());
   SetY(Y() + gamma2*bp*by + gamma*by*T());
   SetZ(Z() + gamma2*bp*bz + gamma*bz*T());
   SetT(gamma*(T() + bp));
}

double Vector4::Rapidity() const
{
   //return rapidity
   return 0.5*log( (E()+Pz()) / (E()-Pz()) );
}

Vector4 &Vector4::operator *= (const LorentzRotation & m)
{
   //multiply this Lorentzvector by m
   return *this = m.VectorMultiplication(*this);
}

Vector4 &Vector4::Transform(const LorentzRotation & m)
{
   //Transform this Lorentzvector
   return *this = m.VectorMultiplication(*this);
}

std::string Vector4::ToString() {
  std::stringstream ss;
  ss << "Vector4(" << fP.X() << "," << fP.Y() << "," << fP.Z() << "," << fE << ")";
  return ss.str();
}

inline double LorentzRotation::XX() const { return fxx; }
inline double LorentzRotation::XY() const { return fxy; }
inline double LorentzRotation::XZ() const { return fxz; }
inline double LorentzRotation::XT() const { return fxt; }
inline double LorentzRotation::YX() const { return fyx; }
inline double LorentzRotation::YY() const { return fyy; }
inline double LorentzRotation::YZ() const { return fyz; }
inline double LorentzRotation::YT() const { return fyt; }
inline double LorentzRotation::ZX() const { return fzx; }
inline double LorentzRotation::ZY() const { return fzy; }
inline double LorentzRotation::ZZ() const { return fzz; }
inline double LorentzRotation::ZT() const { return fzt; }
inline double LorentzRotation::TX() const { return ftx; }
inline double LorentzRotation::TY() const { return fty; }
inline double LorentzRotation::TZ() const { return ftz; }
inline double LorentzRotation::TT() const { return ftt; }

inline LorentzRotation::LorentzRotationRow::LorentzRotationRow
(const LorentzRotation & r, int i) : fRR(&r), fII(i) {}

inline double LorentzRotation::LorentzRotationRow::operator [] (int jj) const {
   return fRR->operator()(fII,jj);
}

inline LorentzRotation::LorentzRotationRow LorentzRotation::operator [] (int i) const {
   return LorentzRotationRow(*this, i);
}

inline LorentzRotation & LorentzRotation::operator = (const LorentzRotation & r) {
   fxx = r.fxx; fxy = r.fxy; fxz = r.fxz; fxt = r.fxt;
   fyx = r.fyx; fyy = r.fyy; fyz = r.fyz; fyt = r.fyt;
   fzx = r.fzx; fzy = r.fzy; fzz = r.fzz; fzt = r.fzt;
   ftx = r.ftx; fty = r.fty; ftz = r.ftz; ftt = r.ftt;
   return *this;
}

//inline LorentzRotation &
//LorentzRotation::operator = (const Rotation & r) {
//  mxx = r.xx(); mxy = r.xy(); mxz = r.xz(); mxt = 0.0;
//  myx = r.yx(); myy = r.yy(); myz = r.yz(); myt = 0.0;
//  mzx = r.zx(); mzy = r.zy(); mzz = r.zz(); mzt = 0.0;
//  mtx = 0.0;    mty = 0.0;    mtz = 0.0;    mtt = 1.0;
//  return *this;
//}

inline LorentzRotation & LorentzRotation::operator = (const Rotation & r) {
   fxx = r.XX(); fxy = r.XY(); fxz = r.XZ(); fxt = 0.0;
   fyx = r.YX(); fyy = r.YY(); fyz = r.YZ(); fyt = 0.0;
   fzx = r.ZX(); fzy = r.ZY(); fzz = r.ZZ(); fzt = 0.0;
   ftx = 0.0;    fty = 0.0;    ftz = 0.0;    ftt = 1.0;
   return *this;
}


//inline bool
//LorentzRotation::operator == (const LorentzRotation & r) const {
//  return (mxx == r.xx() && mxy == r.xy() && mxz == r.xz() && mxt == r.xt() &&
//          myx == r.yx() && myy == r.yy() && myz == r.yz() && myt == r.yt() &&
//          mzx == r.zx() && mzy == r.zy() && mzz == r.zz() && mzt == r.zt() &&
//          mtx == r.tx() && mty == r.ty() && mtz == r.tz() && mtt == r.tt())
//  ? true : false;
//}

inline bool LorentzRotation::operator == (const LorentzRotation & r) const {
   return (fxx == r.fxx && fxy == r.fxy && fxz == r.fxz && fxt == r.fxt &&
           fyx == r.fyx && fyy == r.fyy && fyz == r.fyz && fyt == r.fyt &&
           fzx == r.fzx && fzy == r.fzy && fzz == r.fzz && fzt == r.fzt &&
           ftx == r.ftx && fty == r.fty && ftz == r.ftz && ftt == r.ftt)
   ? true : false;
}

//inline bool
//LorentzRotation::operator != (const LorentzRotation & r) const {
//  return (mxx != r.xx() || mxy != r.xy() || mxz != r.xz() || mxt != r.xt() ||
//          myx != r.yx() || myy != r.yy() || myz != r.yz() || myt != r.yt() ||
//          mzx != r.zx() || mzy != r.zy() || mzz != r.zz() || mzt != r.zt() ||
//          mtx != r.tx() || mty != r.ty() || mtz != r.tz() || mtt != r.tt())
//  ? true : false;
//}

inline bool LorentzRotation::operator != (const LorentzRotation & r) const {
   return (fxx != r.fxx || fxy != r.fxy || fxz != r.fxz || fxt != r.fxt ||
           fyx != r.fyx || fyy != r.fyy || fyz != r.fyz || fyt != r.fyt ||
           fzx != r.fzx || fzy != r.fzy || fzz != r.fzz || fzt != r.fzt ||
           ftx != r.ftx || fty != r.fty || ftz != r.ftz || ftt != r.ftt)
   ? true : false;
}

inline bool LorentzRotation::IsIdentity() const {
   return (fxx == 1.0 && fxy == 0.0 && fxz == 0.0 && fxt == 0.0 &&
           fyx == 0.0 && fyy == 1.0 && fyz == 0.0 && fyt == 0.0 &&
           fzx == 0.0 && fzy == 0.0 && fzz == 1.0 && fzt == 0.0 &&
           ftx == 0.0 && fty == 0.0 && ftz == 0.0 && ftt == 1.0)
   ? true : false;
}


inline Vector4 LorentzRotation::VectorMultiplication(const Vector4 & p) const {
   return Vector4(fxx*p.X()+fxy*p.Y()+fxz*p.Z()+fxt*p.T(),
		  fyx*p.X()+fyy*p.Y()+fyz*p.Z()+fyt*p.T(),
		  fzx*p.X()+fzy*p.Y()+fzz*p.Z()+fzt*p.T(),
		  ftx*p.X()+fty*p.Y()+ftz*p.Z()+ftt*p.T());
}

inline Vector4 LorentzRotation::operator * (const Vector4 & p) const {
   return VectorMultiplication(p);
}

inline LorentzRotation LorentzRotation::operator * (const LorentzRotation & m) const {
   return MatrixMultiplication(m);
}

inline LorentzRotation & LorentzRotation::operator *= (const LorentzRotation & m) {
   return *this = MatrixMultiplication(m);
}

inline LorentzRotation & LorentzRotation::Transform(const LorentzRotation & m) {
   return *this = m.MatrixMultiplication(*this);
}

inline LorentzRotation & LorentzRotation::Transform(const Rotation & m){
   return Transform(LorentzRotation(m));
}

inline LorentzRotation LorentzRotation::Inverse() const {
   return LorentzRotation( fxx,  fyx,  fzx, -ftx,
                            fxy,  fyy,  fzy, -fty,
                            fxz,  fyz,  fzz, -ftz,
                           -fxt, -fyt, -fzt,  ftt);
}

inline LorentzRotation & LorentzRotation::Invert() {
   return *this = Inverse();
}

inline LorentzRotation & LorentzRotation::Boost(double bx, double by, double bz) {
   return Transform(LorentzRotation(bx, by, bz));
}

inline LorentzRotation & LorentzRotation::Boost(const Vector3 & b) {
   return Transform(LorentzRotation(b));
}

inline LorentzRotation & LorentzRotation::RotateX(double angle) {
   return Transform(Rotation().RotateX(angle));
}

inline LorentzRotation & LorentzRotation::RotateY(double angle) {
   return Transform(Rotation().RotateY(angle));
}

inline LorentzRotation & LorentzRotation::RotateZ(double angle) {
   return Transform(Rotation().RotateZ(angle));
}

inline LorentzRotation & LorentzRotation::Rotate(double angle, const Vector3 & axis) {
   return Transform(Rotation().Rotate(angle, axis));
}

inline LorentzRotation & LorentzRotation::Rotate(double angle, const Vector3 * axis) {
   return Transform(Rotation().Rotate(angle, axis));
}


LorentzRotation::LorentzRotation()
  : fxx(1.0), fxy(0.0), fxz(0.0), fxt(0.0),
    fyx(0.0), fyy(1.0), fyz(0.0), fyt(0.0),
    fzx(0.0), fzy(0.0), fzz(1.0), fzt(0.0),
    ftx(0.0), fty(0.0), ftz(0.0), ftt(1.0) {}

LorentzRotation::LorentzRotation(const Rotation & r)
  : fxx(r.XX()), fxy(r.XY()), fxz(r.XZ()), fxt(0.0),
    fyx(r.YX()), fyy(r.YY()), fyz(r.YZ()), fyt(0.0),
    fzx(r.ZX()), fzy(r.ZY()), fzz(r.ZZ()), fzt(0.0),
    ftx(0.0),    fty(0.0),    ftz(0.0),    ftt(1.0) {}

LorentzRotation::LorentzRotation(const LorentzRotation & r) :
    fxx(r.fxx), fxy(r.fxy), fxz(r.fxz), fxt(r.fxt),
    fyx(r.fyx), fyy(r.fyy), fyz(r.fyz), fyt(r.fyt),
    fzx(r.fzx), fzy(r.fzy), fzz(r.fzz), fzt(r.fzt),
    ftx(r.ftx), fty(r.fty), ftz(r.ftz), ftt(r.ftt) {}

LorentzRotation::LorentzRotation(
  double rxx, double rxy, double rxz, double rxt,
  double ryx, double ryy, double ryz, double ryt,
  double rzx, double rzy, double rzz, double rzt,
  double rtx, double rty, double rtz, double rtt)
  : fxx(rxx), fxy(rxy), fxz(rxz), fxt(rxt),
    fyx(ryx), fyy(ryy), fyz(ryz), fyt(ryt),
    fzx(rzx), fzy(rzy), fzz(rzz), fzt(rzt),
    ftx(rtx), fty(rty), ftz(rtz), ftt(rtt) {}

LorentzRotation::LorentzRotation(double bx,
                                   double by,
                                   double bz)
{
   //constructor   
   SetBoost(bx, by, bz);
}

LorentzRotation::LorentzRotation(const Vector3 & p) {
   //copy constructor
   SetBoost(p.X(), p.Y(), p.Z());
}

double LorentzRotation::operator () (int i, int j) const {
   //derefencing operator
   if (i == 0) {
      if (j == 0) { return fxx; }
      if (j == 1) { return fxy; }
      if (j == 2) { return fxz; }
      if (j == 3) { return fxt; }
   } else if (i == 1) {
      if (j == 0) { return fyx; }
      if (j == 1) { return fyy; }
      if (j == 2) { return fyz; }
      if (j == 3) { return fyt; }
   } else if (i == 2) {
      if (j == 0) { return fzx; }
      if (j == 1) { return fzy; }
      if (j == 2) { return fzz; }
      if (j == 3) { return fzt; }
   } else if (i == 3) {
      if (j == 0) { return ftx; }
      if (j == 1) { return fty; }
      if (j == 2) { return ftz; }
      if (j == 3) { return ftt; }
   }
   std::cerr << "LorentzRotation::operator()(i,j): subscripting: bad indices(" << i << "," << j << ")" << std::endl;
   return 0.0;
}

void LorentzRotation::SetBoost(double bx, double by, double bz) {
   //boost this Lorentz vector
   double bp2 = bx*bx + by*by + bz*bz;
   double gamma = 1.0 / sqrt(1.0 - bp2);
   double bgamma = gamma * gamma / (1.0 + gamma);
   fxx = 1.0 + bgamma * bx * bx;
   fyy = 1.0 + bgamma * by * by;
   fzz = 1.0 + bgamma * bz * bz;
   fxy = fyx = bgamma * bx * by;
   fxz = fzx = bgamma * bx * bz;
   fyz = fzy = bgamma * by * bz;
   fxt = ftx = gamma * bx;
   fyt = fty = gamma * by;
   fzt = ftz = gamma * bz;
   ftt = gamma;
}

LorentzRotation LorentzRotation::MatrixMultiplication(const LorentzRotation & b) const {
   //multiply this vector by a matrix
   return LorentzRotation(
    fxx*b.fxx + fxy*b.fyx + fxz*b.fzx + fxt*b.ftx,
    fxx*b.fxy + fxy*b.fyy + fxz*b.fzy + fxt*b.fty,
    fxx*b.fxz + fxy*b.fyz + fxz*b.fzz + fxt*b.ftz,
    fxx*b.fxt + fxy*b.fyt + fxz*b.fzt + fxt*b.ftt,
    fyx*b.fxx + fyy*b.fyx + fyz*b.fzx + fyt*b.ftx,
    fyx*b.fxy + fyy*b.fyy + fyz*b.fzy + fyt*b.fty,
    fyx*b.fxz + fyy*b.fyz + fyz*b.fzz + fyt*b.ftz,
    fyx*b.fxt + fyy*b.fyt + fyz*b.fzt + fyt*b.ftt,
    fzx*b.fxx + fzy*b.fyx + fzz*b.fzx + fzt*b.ftx,
    fzx*b.fxy + fzy*b.fyy + fzz*b.fzy + fzt*b.fty,
    fzx*b.fxz + fzy*b.fyz + fzz*b.fzz + fzt*b.ftz,
    fzx*b.fxt + fzy*b.fyt + fzz*b.fzt + fzt*b.ftt,
    ftx*b.fxx + fty*b.fyx + ftz*b.fzx + ftt*b.ftx,
    ftx*b.fxy + fty*b.fyy + ftz*b.fzy + ftt*b.fty,
    ftx*b.fxz + fty*b.fyz + ftz*b.fzz + ftt*b.ftz,
    ftx*b.fxt + fty*b.fyt + ftz*b.fzt + ftt*b.ftt);
}

std::string LorentzRotation::ToString() {
  std::stringstream ss;
  ss << "LorentzRotation(" << std::setprecision(8) << fxx << " " << fxy << " " << fxz << " " << fxt << std::endl;
  ss << "                " << std::setprecision(8) << fyx << " " << fyy << " " << fyz << " " << fyt << std::endl;
  ss << "                " << std::setprecision(8) << fzx << " " << fzy << " " << fzz << " " << fzt << std::endl;
  ss << "                " << std::setprecision(8) << ftx << " " << fty << " " << ftz << " " << ftt << ")" << std::endl;
  
  return ss.str();
}


#endif
